home *** CD-ROM | disk | FTP | other *** search
/ Tech Arsenal 1 / Tech Arsenal (Arsenal Computer).ISO / tek-19 / gpt32src.zip / CONTOUR.C < prev    next >
C/C++ Source or Header  |  1992-03-25  |  43KB  |  1,244 lines
  1. #ifndef lint
  2. static char *RCSid = "$Id: contour.c,v 3.26 92/03/24 22:35:54 woo Exp Locker: woo $";
  3. #endif
  4.  
  5. /* GNUPLOT - contour.c */
  6. /*
  7.  * Copyright (C) 1986, 1987, 1990, 1991, 1992   Thomas Williams, Colin Kelley
  8.  *
  9.  * Permission to use, copy, and distribute this software and its
  10.  * documentation for any purpose with or without fee is hereby granted, 
  11.  * provided that the above copyright notice appear in all copies and 
  12.  * that both that copyright notice and this permission notice appear 
  13.  * in supporting documentation.
  14.  *
  15.  * Permission to modify the software is granted, but not the right to
  16.  * distribute the modified code.  Modifications are to be distributed 
  17.  * as patches to released version.
  18.  *  
  19.  * This software is provided "as is" without express or implied warranty.
  20.  * 
  21.  *
  22.  * AUTHORS
  23.  * 
  24.  *   Original Software:
  25.  *       Gershon Elber
  26.  * 
  27.  * Send your comments or suggestions to 
  28.  *  info-gnuplot@ames.arc.nasa.gov.
  29.  * This is a mailing list; to join it send a note to 
  30.  *  info-gnuplot-request@ames.arc.nasa.gov.  
  31.  * Send bug reports to
  32.  *  bug-gnuplot@ames.arc.nasa.gov.
  33.  */
  34.  
  35. #include <stdio.h>
  36. #include "plot.h"
  37.  
  38. #define DEFAULT_NUM_OF_ZLEVELS  10  /* Some dflt values (setable via flags). */
  39. #define DEFAULT_NUM_APPROX_PTS  5
  40. #define DEFAULT_BSPLINE_ORDER  3
  41. #define MAX_NUM_OF_ZLEVELS      100 /* Some max. values (setable via flags). */
  42. #define MAX_NUM_APPROX_PTS      100
  43. #define MAX_BSPLINE_ORDER      10
  44.  
  45. #define INTERP_NOTHING   0            /* Kind of interpolations on contours. */
  46. #define INTERP_CUBIC     1                           /* Cubic spline interp. */
  47. #define APPROX_BSPLINE   2                         /* Bspline interpolation. */
  48.  
  49. #define ACTIVE   1                    /* Status of edges at certain Z level. */
  50. #define INACTIVE 2
  51.  
  52. #define OPEN_CONTOUR     1                                 /* Contour kinds. */
  53. #define CLOSED_CONTOUR   2
  54.  
  55. #define EPSILON  1e-5              /* Used to decide if two float are equal. */
  56. #define INFINITY 1e10
  57.  
  58. #ifndef TRUE
  59. #define TRUE     -1
  60. #define FALSE    0
  61. #endif
  62.  
  63. #define DEFAULT_NUM_CONTOURS    10
  64. #define MAX_POINTS_PER_CNTR     100
  65. #define SHIFT_Z_EPSILON        0.000301060 /* Dec. change of poly bndry hit.*/
  66.  
  67. #define abs(x)  ((x) > 0 ? (x) : (-(x)))
  68. #define sqr(x)  ((x) * (x))
  69.  
  70. #ifndef AMIGA_AC_5
  71. extern double sqrt();
  72. #endif /* not AMIGA_AC_5 */
  73. typedef double tri_diag[3];         /* Used to allocate the tri-diag matrix. */
  74. typedef double table_entry[4];           /* Cubic spline interpolation 4 coef. */
  75.  
  76. struct vrtx_struct {
  77.     double X, Y, Z;                       /* The coordinates of this vertex. */
  78.     struct vrtx_struct *next;                             /* To chain lists. */
  79. };
  80.  
  81. struct edge_struct {
  82.     struct poly_struct *poly[2];   /* Each edge belongs to up to 2 polygons. */
  83.     struct vrtx_struct *vertex[2]; /* The two extreme points of this vertex. */
  84.     struct edge_struct *next;                             /* To chain lists. */
  85.     int status, /* Status flag to mark edges in scanning at certain Z level. */
  86.     boundary;                   /* True if this edge is on the boundary. */
  87. };
  88.  
  89. struct poly_struct {
  90.     struct edge_struct *edge[3];           /* As we do triangolation here... */
  91.     struct poly_struct *next;                             /* To chain lists. */
  92. };
  93.  
  94. struct cntr_struct {           /* Contours are saved using this struct list. */
  95.     double X, Y;                          /* The coordinates of this vertex. */
  96.     struct cntr_struct *next;                             /* To chain lists. */
  97. };
  98.  
  99. static int test_boundary;    /* If TRUE look for contours on boundary first. */
  100. static struct gnuplot_contours *contour_list = NULL;
  101. static double crnt_cntr[MAX_POINTS_PER_CNTR * 2];
  102. static int crnt_cntr_pt_index = 0;
  103. static double contour_level = 0.0;
  104. static table_entry *hermit_table = NULL;    /* Hold hermite table constants. */
  105. static int num_of_z_levels = DEFAULT_NUM_OF_ZLEVELS;  /* # Z contour levels. */
  106. static int num_approx_pts = DEFAULT_NUM_APPROX_PTS;/* # pts per approx/inter.*/
  107. static int bspline_order = DEFAULT_BSPLINE_ORDER;   /* Bspline order to use. */
  108. static int interp_kind = INTERP_NOTHING;  /* Linear, Cubic interp., Bspline. */
  109.  
  110. static void gen_contours();
  111. static int update_all_edges();
  112. static struct cntr_struct *gen_one_contour();
  113. static struct cntr_struct *trace_contour();
  114. static struct cntr_struct *update_cntr_pt();
  115. static int fuzzy_equal();
  116. static void gen_triangle();
  117. static struct vrtx_struct *gen_vertices();
  118. static struct edge_struct *gen_edges_middle();
  119. static struct edge_struct *gen_edges();
  120. static struct poly_struct *gen_polys();
  121. static void free_contour();
  122. static void put_contour();
  123. static put_contour_nothing();
  124. static put_contour_cubic();
  125. static put_contour_bspline();
  126. static calc_tangent();
  127. static int count_contour();
  128. static complete_spline_interp();
  129. static calc_hermit_table();
  130. static hermit_interp();
  131. static prepare_spline_interp();
  132. static int solve_tri_diag();
  133. static gen_bspline_approx();
  134. static double fetch_knot();
  135. static eval_bspline();
  136.  
  137. /*
  138.  * Entry routine to this whole set of contouring module.
  139.  */
  140. struct gnuplot_contours *contour(num_isolines, iso_lines,
  141.                  ZLevels, approx_pts, kind, order1)
  142. int num_isolines;
  143. struct iso_curve *iso_lines;
  144. int ZLevels, approx_pts, kind, order1;
  145. {
  146.     int i;
  147.     struct poly_struct *p_polys, *p_poly;
  148.     struct edge_struct *p_edges, *p_edge;
  149.     struct vrtx_struct *p_vrts, *p_vrtx;
  150.     double x_min, y_min, z_min, x_max, y_max, z_max, z, dz, z_scale = 1.0;
  151.  
  152.     num_of_z_levels = ZLevels;
  153.     num_approx_pts = approx_pts;
  154.     bspline_order = order1 - 1;
  155.     interp_kind = kind;
  156.  
  157.     contour_list = NULL;
  158.  
  159.     if (interp_kind == INTERP_CUBIC) calc_hermit_table();
  160.  
  161.     gen_triangle(num_isolines, iso_lines, &p_polys, &p_edges, &p_vrts,
  162.         &x_min, &y_min, &z_min, &x_max, &y_max, &z_max);
  163.     dz = (z_max - z_min) / (num_of_z_levels+1);
  164.     /* Step from z_min+dz upto z_max-dz in num_of_z_levels times. */
  165.     z = z_min + dz;
  166.     crnt_cntr_pt_index = 0;
  167.  
  168.     for (i=0; i<num_of_z_levels; i++, z += dz) {
  169.     contour_level = z;
  170.     gen_contours(p_edges, z + dz * SHIFT_Z_EPSILON, x_min, x_max,
  171.                             y_min, y_max);
  172.     }
  173.  
  174.     /* Free all contouring related temporary data. */
  175.     while (p_polys) {
  176.     p_poly = p_polys -> next;
  177.     free (p_polys);
  178.     p_polys = p_poly;
  179.     }
  180.     while (p_edges) {
  181.     p_edge = p_edges -> next;
  182.     free (p_edges);
  183.     p_edges = p_edge;
  184.     }
  185.     while (p_vrts) {
  186.     p_vrtx = p_vrts -> next;
  187.     free (p_vrts);
  188.     p_vrts = p_vrtx;
  189.     }
  190.  
  191.     if (interp_kind == INTERP_CUBIC) free(hermit_table);
  192.  
  193.     return contour_list;
  194. }
  195.  
  196. /*
  197.  * Adds another point to the currently build contour.
  198.  */
  199. add_cntr_point(x, y)
  200. double x, y;
  201. {
  202.     int index;
  203.  
  204.     if (crnt_cntr_pt_index >= MAX_POINTS_PER_CNTR-1) {
  205.     index = crnt_cntr_pt_index - 1;
  206.     end_crnt_cntr();
  207.     crnt_cntr[0] = crnt_cntr[index * 2];
  208.     crnt_cntr[1] = crnt_cntr[index * 2 + 1];
  209.     crnt_cntr_pt_index = 1; /* Keep the last point as first of this one. */
  210.     }
  211.     crnt_cntr[crnt_cntr_pt_index * 2] = x;
  212.     crnt_cntr[crnt_cntr_pt_index * 2 + 1] = y;
  213.     crnt_cntr_pt_index++;
  214. }
  215.  
  216. /*
  217.  * Done with current contour - create gnuplot data structure for it.
  218.  */
  219. end_crnt_cntr()
  220. {
  221.     int i;
  222.     struct gnuplot_contours *cntr = (struct gnuplot_contours *)
  223.                     alloc(sizeof(struct gnuplot_contours),
  224.                           "gnuplot_contour");
  225.  
  226.     cntr->coords = (struct coordinate *) alloc(sizeof(struct coordinate) *
  227.                               crnt_cntr_pt_index,
  228.                            "contour coords");
  229.     for (i=0; i<crnt_cntr_pt_index; i++) {
  230.     cntr->coords[i].x = crnt_cntr[i * 2];
  231.     cntr->coords[i].y = crnt_cntr[i * 2 + 1];
  232.     cntr->coords[i].z = contour_level;
  233.     }
  234.     cntr->num_pts = crnt_cntr_pt_index;
  235.  
  236.     cntr->next = contour_list;
  237.     contour_list = cntr;
  238.  
  239.     crnt_cntr_pt_index = 0;
  240. }
  241.  
  242. /*
  243.  * Generates all contours by tracing the intersecting triangles.
  244.  */
  245. static void gen_contours(p_edges, z_level, x_min, x_max, y_min, y_max)
  246. struct edge_struct *p_edges;
  247. double z_level, x_min, x_max, y_min, y_max;
  248. {
  249.     int num_active,                        /* Number of edges marked ACTIVE. */
  250.     contour_kind;                /* One of OPEN_CONTOUR, CLOSED_CONTOUR. */
  251.     struct cntr_struct *p_cntr;
  252.  
  253.     num_active = update_all_edges(p_edges, z_level);           /* Do pass 1. */
  254.  
  255.     test_boundary = TRUE;        /* Start to look for contour on boundaries. */
  256.  
  257.     while (num_active > 0) {                                   /* Do Pass 2. */
  258.         /* Generate One contour (and update MumActive as needed): */
  259.     p_cntr = gen_one_contour(p_edges, z_level, &contour_kind, &num_active);
  260.     put_contour(p_cntr, z_level, x_min, x_max, y_min, y_max,
  261.                   contour_kind); /* Emit it in requested format. */
  262.     }
  263. }
  264.  
  265. /*
  266.  * Does pass 1, or marks the edges which are active (crosses this z_level)
  267.  * as ACTIVE, and the others as INACTIVE:
  268.  * Returns number of active edges (marked ACTIVE).
  269.  */
  270. static int update_all_edges(p_edges, z_level)
  271. struct edge_struct *p_edges;
  272. double z_level;
  273. {
  274.     int count = 0;
  275.  
  276.     while (p_edges) {
  277.     if (((p_edges -> vertex[0] -> Z >= z_level) &&
  278.          (p_edges -> vertex[1] -> Z <= z_level)) ||
  279.         ((p_edges -> vertex[1] -> Z >= z_level) &&
  280.          (p_edges -> vertex[0] -> Z <= z_level))) {
  281.         p_edges -> status = ACTIVE;
  282.         count++;
  283.     }
  284.     else p_edges -> status = INACTIVE;
  285.     p_edges = p_edges -> next;
  286.     }
  287.  
  288.     return count;
  289. }
  290.  
  291. /*
  292.  * Does pass 2, or find one complete contour out of the triangolation data base:
  293.  * Returns a pointer to the contour (as linked list), contour_kind is set to
  294.  * one of OPEN_CONTOUR, CLOSED_CONTOUR, and num_active is updated.
  295.  */
  296. static struct cntr_struct *gen_one_contour(p_edges, z_level, contour_kind,
  297.                                 num_active)
  298. struct edge_struct *p_edges;
  299. double z_level;
  300. int *contour_kind, *num_active;
  301. {
  302.     struct edge_struct *pe_temp;
  303.  
  304.     if (test_boundary) {    /* Look for something to start with on boundary: */
  305.     pe_temp = p_edges;
  306.     while (pe_temp) {
  307.         if ((pe_temp -> status == ACTIVE) && (pe_temp -> boundary)) break;
  308.         pe_temp = pe_temp -> next;
  309.     }
  310.     if (!pe_temp) test_boundary = FALSE;/* No more contours on boundary. */
  311.     else {
  312.         *contour_kind = OPEN_CONTOUR;
  313.         return trace_contour(pe_temp, z_level, num_active, *contour_kind);
  314.     }
  315.     }
  316.  
  317.     if (!test_boundary) {        /* Look for something to start with inside: */
  318.     pe_temp = p_edges;
  319.     while (pe_temp) {
  320.         if ((pe_temp -> status == ACTIVE) && (!(pe_temp -> boundary)))
  321.         break;
  322.         pe_temp = pe_temp -> next;
  323.     }
  324.     if (!pe_temp) {
  325.         *num_active = 0;
  326.         return NULL;
  327.     }
  328.     else {
  329.         *contour_kind = CLOSED_CONTOUR;
  330.         return trace_contour(pe_temp, z_level, num_active, *contour_kind);
  331.     }
  332.     }
  333.     return NULL;             /* We should never be here, but lint... */
  334. }
  335.  
  336. /*
  337.  * Search the data base along a contour starts at the edge pe_start until
  338.  * a boundary edge is detected or until we close the loop back to pe_start.
  339.  * Returns a linked list of all the points on the contour
  340.  * Also decreases num_active by the number of points on contour.
  341.  */
  342. static struct cntr_struct *trace_contour(pe_start, z_level, num_active,
  343.                                 contour_kind)
  344. struct edge_struct *pe_start;
  345. double z_level;
  346. int *num_active, contour_kind;
  347. {
  348.     int i, in_middle;       /* If TRUE the z_level is in the middle of edge. */
  349.     struct cntr_struct *p_cntr, *pc_tail;
  350.     struct edge_struct *p_edge = pe_start, *p_next_edge;
  351.     struct poly_struct *p_poly, *PLastpoly = NULL;
  352.  
  353.     /* Generate the header of the contour - the point on pe_start. */
  354.     if (contour_kind == OPEN_CONTOUR) pe_start -> status = INACTIVE;
  355.     (*num_active)--;
  356.     p_cntr = pc_tail = update_cntr_pt(pe_start, z_level, &in_middle);
  357.     if (!in_middle) {
  358.     return NULL;
  359.     }
  360.  
  361.     do {
  362.     /* Find polygon to continue (Not where we came from - PLastpoly): */
  363.     if (p_edge -> poly[0] == PLastpoly) p_poly = p_edge -> poly[1];
  364.     else p_poly = p_edge -> poly[0];
  365.     p_next_edge = NULL;          /* In case of error, remains NULL. */
  366.     for (i=0; i<3; i++)              /* Test the 3 edges of the polygon: */
  367.         if (p_poly -> edge[i] != p_edge)
  368.             if (p_poly -> edge[i] -> status == ACTIVE)
  369.             p_next_edge = p_poly -> edge[i];
  370.     if (!p_next_edge) {
  371.         pc_tail -> next = NULL;
  372.         free_contour(p_cntr);
  373.         return NULL;
  374.     }
  375.     p_edge = p_next_edge;
  376.     PLastpoly = p_poly;
  377.     p_edge -> status = INACTIVE;
  378.     (*num_active)--;
  379.     pc_tail -> next = update_cntr_pt(p_edge, z_level, &in_middle);
  380.     if (!in_middle) {
  381.         pc_tail -> next = NULL;
  382.         free_contour(p_cntr);
  383.         return NULL;
  384.     }
  385.         pc_tail = pc_tail -> next;
  386.     }
  387.     while ((pe_start != p_edge) && (!p_edge -> boundary));
  388.     pc_tail -> next = NULL;
  389.  
  390.     return p_cntr;
  391. }
  392.  
  393. /*
  394.  * Allocates one contour location and update it to to correct position
  395.  * according to z_level and edge p_edge. if z_level is found to be at
  396.  * one of the extreme points nothing is allocated (NULL is returned)
  397.  * and in_middle is set to FALSE.
  398.  */
  399. static struct cntr_struct *update_cntr_pt(p_edge, z_level, in_middle)
  400. struct edge_struct *p_edge;
  401. double z_level;
  402. int *in_middle;
  403. {
  404.     double t;
  405.     struct cntr_struct *p_cntr;
  406.  
  407.     t = (z_level - p_edge -> vertex[0] -> Z) /
  408.     (p_edge -> vertex[1] -> Z - p_edge -> vertex[0] -> Z);
  409.  
  410.     if (fuzzy_equal(t, 1.0) || fuzzy_equal(t, 0.0)) {
  411.         *in_middle = FALSE;
  412.         return NULL;
  413.     }
  414.     else {
  415.     *in_middle = TRUE;
  416.     p_cntr = (struct cntr_struct *) alloc(sizeof(struct cntr_struct),
  417.                             "contour cntr_struct");
  418.     p_cntr -> X = p_edge -> vertex[1] -> X * t +
  419.              p_edge -> vertex[0] -> X * (1-t);
  420.     p_cntr -> Y = p_edge -> vertex[1] -> Y * t +
  421.              p_edge -> vertex[0] -> Y * (1-t);
  422.     return p_cntr;
  423.     }
  424. }
  425.  
  426. /*
  427.  * Simple routine to decide if two real values are equal by simply
  428.  * calculating the relative/absolute error between them (< EPSILON).
  429.  */
  430. static int fuzzy_equal(x, y)
  431. double x, y;
  432. {
  433.     if (abs(x) > EPSILON)            /* Calculate relative error: */
  434.         return (abs((x - y) / x) < EPSILON);
  435.     else                    /* Calculate absolute error: */
  436.     return (abs(x - y) < EPSILON);
  437. }
  438.  
  439. /*
  440.  * Generate the triangles.
  441.  * Returns the lists (vrtxs edges & polys) via pointers to their heads.
  442.  */
  443. static void gen_triangle(num_isolines, iso_lines, p_polys, p_edges,
  444.     p_vrts, x_min, y_min, z_min, x_max, y_max, z_max)
  445. int num_isolines;
  446. struct iso_curve *iso_lines;
  447. struct poly_struct **p_polys;
  448. struct edge_struct **p_edges;
  449. struct vrtx_struct **p_vrts;
  450. double *x_min, *y_min, *z_min, *x_max, *y_max, *z_max;
  451. {
  452.     int i, grid_x_max = iso_lines->p_count;
  453.     struct vrtx_struct *p_vrtx1, *p_vrtx2, *pv_temp;
  454.     struct edge_struct *p_edge1, *p_edge2, *pe_tail1, *pe_tail2, *pe_temp,
  455.               *p_edge_middle, *pe_m_tail;
  456.     struct poly_struct *p_poly, *pp_tail;
  457.  
  458.     *p_polys = NULL;
  459.     *p_edges = NULL;
  460.     *p_vrts = NULL;
  461.     *z_min = INFINITY;
  462.     *y_min = INFINITY;
  463.     *x_min = INFINITY;
  464.     *z_max = -INFINITY;
  465.     *y_max = -INFINITY;
  466.     *x_max = -INFINITY;
  467.  
  468.     /* Read 1st row. */
  469.     p_vrtx1 = gen_vertices(grid_x_max, iso_lines->points,
  470.                x_min, y_min, z_min, x_max, y_max, z_max);
  471.     *p_vrts = p_vrtx1;
  472.     /* Gen. its edges.*/
  473.     pe_temp = p_edge1 = gen_edges(grid_x_max, p_vrtx1, &pe_tail1);
  474.     for (i = 1; i < grid_x_max; i++) {/* Mark one side of edges as boundary. */
  475.     pe_temp -> poly[1] = NULL;
  476.     pe_temp = pe_temp -> next;
  477.     }
  478.     for (i = 1; i < num_isolines; i++) { /* Read next column and gen. polys. */
  479.     iso_lines = iso_lines->next;
  480.         /* Get row into list. */
  481.         p_vrtx2 = gen_vertices(grid_x_max, iso_lines->points,
  482.                    x_min, y_min, z_min, x_max, y_max, z_max);
  483.         /* Generate its edges. */
  484.         p_edge2 = gen_edges(grid_x_max, p_vrtx2, &pe_tail2);
  485.     /* Generate edges from one vertex list to the other one: */
  486.     p_edge_middle = gen_edges_middle(grid_x_max, p_vrtx1, p_vrtx2,
  487.                                  &pe_m_tail);
  488.  
  489.     /* Now we can generate the polygons themselves (triangles). */
  490.     p_poly = gen_polys(grid_x_max, p_edge1, p_edge_middle, p_edge2,
  491.                                  &pp_tail);
  492.         pe_tail1 -> next = (*p_edges);      /* Chain new edges to main list. */
  493.         pe_m_tail -> next = p_edge1;
  494.     *p_edges = p_edge_middle;
  495.     pe_tail1 = pe_tail2;
  496.     p_edge1 = p_edge2;
  497.  
  498.     pv_temp = p_vrtx2;
  499.     while (pv_temp -> next) pv_temp = pv_temp -> next;
  500.     pv_temp -> next = *p_vrts;
  501.     *p_vrts = p_vrtx1 = p_vrtx2;
  502.  
  503.         pp_tail -> next = (*p_polys);       /* Chain new polys to main list. */
  504.     *p_polys = p_poly;
  505.     }
  506.  
  507.     pe_temp = p_edge1;
  508.     for (i = 1; i < grid_x_max; i++) {/* Mark one side of edges as boundary. */
  509.     pe_temp -> poly[0] = NULL;
  510.     pe_temp = pe_temp -> next;
  511.     }
  512.  
  513.     pe_tail1 -> next = (*p_edges);    /* Chain last edges list to main list. */
  514.     *p_edges = p_edge1;
  515.  
  516.     /* Update the boundary flag, saved in each edge, and update indexes: */
  517.     pe_temp = (*p_edges);
  518.     i = 1;
  519.  
  520.     while (pe_temp) {
  521.     pe_temp -> boundary = (!(pe_temp -> poly[0])) ||
  522.                   (!(pe_temp -> poly[1]));
  523.     pe_temp = pe_temp -> next;
  524.     }
  525. }
  526.  
  527. /*
  528.  * Handles grid_x_max 3D points (One row) and generate linked list for them.
  529.  */
  530. static struct vrtx_struct *gen_vertices(grid_x_max, points,
  531.                       x_min, y_min, z_min, x_max, y_max, z_max)
  532. int grid_x_max;
  533. struct coordinate *points;
  534. double *x_min, *y_min, *z_min, *x_max, *y_max, *z_max;
  535. {
  536.     int i;
  537.     struct vrtx_struct *p_vrtx, *pv_tail, *pv_temp;
  538.  
  539.     for (i=0; i<grid_x_max; i++) {/* Get a point and generate the structure. */
  540.         pv_temp = (struct vrtx_struct *) alloc(sizeof(struct vrtx_struct),
  541.                         "contour vertex");
  542.     pv_temp -> X = points[i].x;
  543.     pv_temp -> Y = points[i].y;
  544.     pv_temp -> Z = points[i].z;
  545.  
  546.     if (pv_temp -> X > *x_max) *x_max = pv_temp -> X; /* Update min/max. */
  547.     if (pv_temp -> Y > *y_max) *y_max = pv_temp -> Y;
  548.     if (pv_temp -> Z > *z_max) *z_max = pv_temp -> Z;
  549.     if (pv_temp -> X < *x_min) *x_min = pv_temp -> X;
  550.     if (pv_temp -> Y < *y_min) *y_min = pv_temp -> Y;
  551.     if (pv_temp -> Z < *z_min) *z_min = pv_temp -> Z;
  552.  
  553.     if (i == 0)                              /* First vertex in row: */
  554.         p_vrtx = pv_tail = pv_temp;
  555.     else {
  556.         pv_tail -> next = pv_temp;   /* Stick new record as last one. */
  557.         pv_tail = pv_tail -> next;    /* And continue to last record. */
  558.     }
  559.     }
  560.     pv_tail -> next = NULL;
  561.  
  562.     return p_vrtx;
  563. }
  564.  
  565. /*
  566.  * Combines N vertices in pair to form N-1 edges.
  567.  * Returns pointer to the edge list (pe_tail will point on last edge in list).
  568.  */
  569. static struct edge_struct *gen_edges(grid_x_max, p_vrtx, pe_tail)
  570. int grid_x_max;
  571. struct vrtx_struct *p_vrtx;
  572. struct edge_struct **pe_tail;
  573. {
  574.     int i;
  575.     struct edge_struct *p_edge, *pe_temp;
  576.  
  577.     for (i=0; i<grid_x_max-1; i++) {         /* Generate grid_x_max-1 edges: */
  578.     pe_temp = (struct edge_struct *) alloc(sizeof(struct edge_struct),
  579.                         "contour edge");
  580.     pe_temp -> vertex[0] = p_vrtx;              /* First vertex of edge. */
  581.     p_vrtx = p_vrtx -> next;                     /* Skip to next vertex. */
  582.     pe_temp -> vertex[1] = p_vrtx;             /* Second vertex of edge. */
  583.         if (i == 0)                                    /* First edge in row: */
  584.         p_edge = (*pe_tail) = pe_temp;
  585.     else {
  586.         (*pe_tail) -> next = pe_temp;   /* Stick new record as last one. */
  587.         *pe_tail = (*pe_tail) -> next;   /* And continue to last record. */
  588.      }
  589.     }
  590.     (*pe_tail) -> next = NULL;
  591.  
  592.     return p_edge;
  593. }
  594.  
  595. /*
  596.  * Combines 2 lists of N vertices each into edge list:
  597.  * The dots (.) are the vertices list, and the              .  .  .  .
  598.  *  edges generated are alternations of vertical edges      |\ |\ |\ |
  599.  *  (|) and diagonal ones (\).                              | \| \| \|
  600.  *  A pointer to edge list (alternate | , \) is returned    .  .  .  .
  601.  * Note this list will have (2*grid_x_max-1) edges (pe_tail points on last
  602.  * record).
  603.  */
  604. static struct edge_struct *gen_edges_middle(grid_x_max, p_vrtx1, p_vrtx2,
  605.                                 pe_tail)
  606. int grid_x_max;
  607. struct vrtx_struct *p_vrtx1, *p_vrtx2;
  608. struct edge_struct **pe_tail;
  609. {
  610.     int i;
  611.     struct edge_struct *p_edge, *pe_temp;
  612.  
  613.     /* Gen first (|). */
  614.     pe_temp = (struct edge_struct *) alloc(sizeof(struct edge_struct),
  615.                             "contour edge");
  616.     pe_temp -> vertex[0] = p_vrtx2;                 /* First vertex of edge. */
  617.     pe_temp -> vertex[1] = p_vrtx1;                /* Second vertex of edge. */
  618.     p_edge = (*pe_tail) = pe_temp;
  619.  
  620.     /* Advance in vrtx list grid_x_max-1 times, and gen. 2 edges /| for each.*/
  621.     for (i=0; i<grid_x_max-1; i++) {
  622.     /* The / edge. */
  623.     pe_temp = (struct edge_struct *) alloc(sizeof(struct edge_struct),
  624.                             "contour edge");
  625.     pe_temp -> vertex[0] = p_vrtx1;             /* First vertex of edge. */
  626.     pe_temp -> vertex[1] = p_vrtx2 -> next;    /* Second vertex of edge. */
  627.         (*pe_tail) -> next = pe_temp;       /* Stick new record as last one. */
  628.     *pe_tail = (*pe_tail) -> next;       /* And continue to last record. */
  629.  
  630.     /* The | edge. */
  631.     pe_temp = (struct edge_struct *) alloc(sizeof(struct edge_struct),
  632.                             "contour edge");
  633.     pe_temp -> vertex[0] = p_vrtx2 -> next;     /* First vertex of edge. */
  634.     pe_temp -> vertex[1] = p_vrtx1 -> next;    /* Second vertex of edge. */
  635.         (*pe_tail) -> next = pe_temp;       /* Stick new record as last one. */
  636.     *pe_tail = (*pe_tail) -> next;       /* And continue to last record. */
  637.  
  638.         p_vrtx1 = p_vrtx1 -> next;   /* Skip to next vertices in both lists. */
  639.         p_vrtx2 = p_vrtx2 -> next;
  640.     }
  641.     (*pe_tail) -> next = NULL;
  642.  
  643.     return p_edge;
  644. }
  645.  
  646. /*
  647.  * Combines 3 lists of edges into triangles:
  648.  * 1. p_edge1: Top horizontal edge list:        -----------------------
  649.  * 2. p_edge_middge: middle edge list:         |\  |\  |\  |\  |\  |\  |
  650.  *                                             |  \|  \|  \|  \|  \|  \|
  651.  * 3. p_edge2: Bottom horizontal edge list:     -----------------------
  652.  * Note that p_edge1/2 lists has grid_x_max-1 edges, while p_edge_middle has
  653.  * (2*grid_x_max-1) edges.
  654.  * The routine simple scans the two list    Upper 1         Lower
  655.  * and generate two triangle upper one        ----         | \
  656.  * and lower one from the lists:             0\   |2      0|   \1
  657.  * (Nums. are edges order in polys)             \ |         ----
  658.  * The routine returns a pointer to a                         2
  659.  * polygon list (pp_tail points on last polygon).          1
  660.  *                                                   -----------
  661.  * In addition, the edge lists are updated -        | \   0     |
  662.  * each edge has two pointers on the two            |   \       |
  663.  * (one active if boundary) polygons which         0|1   0\1   0|1
  664.  * uses it. These two pointer to polygons           |       \   |
  665.  * are named: poly[0], poly[1]. The diagram         |    1    \ |
  666.  * on the right show how they are used for the       -----------
  667.  * upper and lower polygons.                             0
  668.  */
  669. static struct poly_struct *gen_polys(grid_x_max, p_edge1, p_edge_middle,
  670.                             p_edge2, pp_tail)
  671. int grid_x_max;
  672. struct edge_struct *p_edge1, *p_edge_middle, *p_edge2;
  673. struct poly_struct **pp_tail;
  674. {
  675.     int i;
  676.     struct poly_struct *p_poly, *pp_temp;
  677.  
  678.     p_edge_middle -> poly[0] = NULL;                /* Its boundary! */
  679.  
  680.     /* Advance in vrtx list grid_x_max-1 times, and gen. 2 polys for each. */
  681.     for (i=0; i<grid_x_max-1; i++) {
  682.     /* The Upper. */
  683.     pp_temp = (struct poly_struct *) alloc(sizeof(struct poly_struct),
  684.                             "contour poly");
  685.     /* Now update polys about its edges, and edges about the polygon. */
  686.     pp_temp -> edge[0] = p_edge_middle -> next;
  687.     p_edge_middle -> next -> poly[1] = pp_temp;
  688.     pp_temp -> edge[1] = p_edge1;
  689.     p_edge1 -> poly[0] = pp_temp;
  690.     pp_temp -> edge[2] = p_edge_middle -> next -> next;
  691.     p_edge_middle -> next -> next -> poly[0] = pp_temp;
  692.     if (i == 0)                   /* Its first one in list: */
  693.         p_poly = (*pp_tail) = pp_temp;
  694.     else {
  695.         (*pp_tail) -> next = pp_temp;
  696.         *pp_tail = (*pp_tail) -> next;
  697.     }
  698.  
  699.     /* The Lower. */
  700.     pp_temp = (struct poly_struct *) alloc(sizeof(struct poly_struct),
  701.                             "contour poly");
  702.     /* Now update polys about its edges, and edges about the polygon. */
  703.     pp_temp -> edge[0] = p_edge_middle;
  704.     p_edge_middle -> poly[1] = pp_temp;
  705.     pp_temp -> edge[1] = p_edge_middle -> next;
  706.     p_edge_middle -> next -> poly[0] = pp_temp;
  707.     pp_temp -> edge[2] = p_edge2;
  708.     p_edge2 -> poly[1] = pp_temp;
  709.     (*pp_tail) -> next = pp_temp;
  710.     *pp_tail = (*pp_tail) -> next;
  711.  
  712.         p_edge1 = p_edge1 -> next;
  713.         p_edge2 = p_edge2 -> next;
  714.         p_edge_middle = p_edge_middle -> next -> next;
  715.     }
  716.     p_edge_middle -> poly[1] = NULL;                /* Its boundary! */
  717.     (*pp_tail) -> next = NULL;
  718.  
  719.     return p_poly;
  720. }
  721.  
  722. /*
  723.  * Calls the (hopefully) desired interpolation/approximation routine.
  724.  */
  725. static void put_contour(p_cntr, z_level, x_min, x_max, y_min, y_max, contr_kind)
  726. struct cntr_struct *p_cntr;
  727. double z_level, x_min, x_max, y_min, y_max;
  728. int contr_kind;
  729. {
  730.     if (!p_cntr) return;            /* Nothing to do if it is empty contour. */
  731.  
  732.     switch (interp_kind) {
  733.     case INTERP_NOTHING:              /* No interpolation/approximation. */
  734.         put_contour_nothing(p_cntr);
  735.         break;
  736.     case INTERP_CUBIC:                    /* Cubic spline interpolation. */
  737.         put_contour_cubic(p_cntr, z_level, x_min, x_max, y_min, y_max,
  738.                                 contr_kind);
  739.         break;
  740.     case APPROX_BSPLINE:                       /* Bspline approximation. */
  741.         put_contour_bspline(p_cntr, z_level, x_min, x_max, y_min, y_max,
  742.                                     contr_kind);
  743.         break;
  744.     }
  745.  
  746.     free_contour(p_cntr);
  747. }
  748.  
  749. /*
  750.  * Simply puts contour coordinates in order with no interpolation or
  751.  * approximation.
  752.  */
  753. static put_contour_nothing(p_cntr)
  754. struct cntr_struct *p_cntr;
  755. {
  756.     while (p_cntr) {
  757.     add_cntr_point(p_cntr -> X, p_cntr -> Y);
  758.     p_cntr = p_cntr -> next;
  759.     }
  760.     end_crnt_cntr();
  761. }
  762.  
  763. /*
  764.  * Find Complete Cubic Spline Interpolation.
  765.  */
  766. static put_contour_cubic(p_cntr, z_level, x_min, x_max, y_min, y_max,
  767.                                  contr_kind)
  768. struct cntr_struct *p_cntr;
  769. double z_level, x_min, x_max, y_min, y_max;
  770. int contr_kind;
  771. {
  772.     int num_pts, i;
  773.     double tx1, ty1, tx2, ty2;                    /* Tangents at end points. */
  774.     struct cntr_struct *pc_temp;
  775.  
  776.     num_pts = count_contour(p_cntr);         /* Number of points in contour. */
  777.  
  778.     if (num_pts > 2) {  /* Take into account 3 points in tangent estimation. */
  779.     calc_tangent(3, p_cntr -> X, p_cntr -> next -> X,
  780.             p_cntr -> next -> next -> X,
  781.             p_cntr -> Y, p_cntr -> next -> Y,
  782.             p_cntr -> next -> next -> Y, &tx1, &ty1);
  783.     pc_temp = p_cntr;
  784.     for (i=3; i<num_pts; i++) pc_temp = pc_temp -> next;/* Go to the end.*/
  785.     calc_tangent(3, pc_temp -> next -> next -> X,
  786.              pc_temp -> next -> X, pc_temp -> X,
  787.             pc_temp -> next -> next -> Y,
  788.             pc_temp -> next -> Y, pc_temp -> Y, &tx2, &ty2);
  789.         tx2 = (-tx2);   /* Inverse the vector as we need opposite direction. */
  790.         ty2 = (-ty2);
  791.     }
  792.     /* If following (num_pts > 1) is TRUE then exactly 2 points in contour.  */
  793.     else if (num_pts > 1) {/* Take into account 2 points in tangent estimat. */
  794.     calc_tangent(2, p_cntr -> X, p_cntr -> next -> X, 0.0,
  795.             p_cntr -> Y, p_cntr -> next -> Y, 0.0, &tx1, &ty1);
  796.     calc_tangent(2, p_cntr -> next -> X, p_cntr -> X, 0.0,
  797.             p_cntr -> next -> Y, p_cntr -> Y, 0.0, &tx2, &ty2);
  798.         tx2 = (-tx2);   /* Inverse the vector as we need opposite direction. */
  799.         ty2 = (-ty2);
  800.     }
  801.     else return;            /* Only one point (???) - ignore it. */
  802.  
  803.     switch (contr_kind) {
  804.     case OPEN_CONTOUR:
  805.         break;
  806.     case CLOSED_CONTOUR:
  807.         tx1 = tx2 = (tx1 + tx2) / 2.0;         /* Make tangents equal. */
  808.         ty1 = ty2 = (ty1 + ty2) / 2.0;
  809.         break;
  810.     }
  811.     complete_spline_interp(p_cntr, num_pts, 0.0, 1.0, tx1, ty1, tx2, ty2);
  812.     end_crnt_cntr();
  813. }
  814.  
  815. /*
  816.  * Find Bspline approximation for this data set.
  817.  * Uses global variable num_approx_pts to determine number of samples per
  818.  * interval, where the knot vector intervals are assumed to be uniform, and
  819.  * Global variable bspline_order for the order of Bspline to use.
  820.  */
  821. static put_contour_bspline(p_cntr, z_level, x_min, x_max, y_min, y_max,
  822.                                 contr_kind)
  823. struct cntr_struct *p_cntr;
  824. double z_level, x_min, x_max,  y_min, y_max;
  825. int contr_kind;
  826. {
  827.     int num_pts, i, order = bspline_order;
  828.     struct cntr_struct *pc_temp;
  829.  
  830.     num_pts = count_contour(p_cntr);         /* Number of points in contour. */
  831.     if (num_pts < 2) return;     /* Cannt do nothing if empty or one points! */
  832.     /* Order must be less than number of points in curve - fix it if needed. */
  833.     if (order > num_pts - 1) order = num_pts - 1;
  834.  
  835.     gen_bspline_approx(p_cntr, num_pts, order, contr_kind);
  836.     end_crnt_cntr();
  837. }
  838.  
  839. /*
  840.  * Estimate the tangents according to the n last points where n might be
  841.  * 2 or 3 (if 2 onlt x1, x2).
  842.  */
  843. static calc_tangent(n, x1, x2, x3, y1, y2, y3, tx, ty)
  844. int n;
  845. double x1, x2, x3, y1, y2, y3, *tx, *ty;
  846. {
  847.     double v1[2], v2[2], v1_magnitude, v2_magnitude;
  848.  
  849.     switch (n) {
  850.     case 2:
  851.         *tx = (x2 - x1) * 0.3;
  852.         *ty = (y2 - y1) * 0.3;
  853.         break;
  854.     case 3:
  855.         v1[0] = x2 - x1;   v1[1] = y2 - y1;
  856.         v2[0] = x3 - x2;   v2[1] = y3 - y2;
  857.         v1_magnitude = sqrt(sqr(v1[0]) + sqr(v1[1]));
  858.         v2_magnitude = sqrt(sqr(v2[0]) + sqr(v2[1]));
  859.         *tx = (v1[0] / v1_magnitude) - (v2[0] / v2_magnitude) * 0.1;
  860.         *tx *= v1_magnitude * 0.1;  /* Make tangent less than magnitude. */
  861.         *ty = (v1[1] / v1_magnitude) - (v2[1] / v2_magnitude) * 0.1;
  862.         *ty *= v1_magnitude * 0.1;  /* Make tangent less than magnitude. */
  863.         break;
  864.     default:                       /* Should not happen! */
  865.         (*ty) = 0.1;
  866.         *tx = 0.1;
  867.         break;
  868.     }
  869. }
  870.  
  871. /*
  872.  * Free all elements in the contour list.
  873.  */
  874. static void free_contour(p_cntr)
  875. struct cntr_struct *p_cntr;
  876. {
  877.     struct cntr_struct *pc_temp;
  878.  
  879.     while (p_cntr) {
  880.     pc_temp = p_cntr;
  881.     p_cntr = p_cntr -> next;
  882.     free((char *) pc_temp);
  883.     }
  884. }
  885.  
  886. /*
  887.  * Counts number of points in contour.
  888.  */
  889. static int count_contour(p_cntr)
  890. struct cntr_struct *p_cntr;
  891. {
  892.     int count = 0;
  893.  
  894.     while (p_cntr) {
  895.     count++;
  896.     p_cntr = p_cntr -> next;
  897.     }
  898.     return count;
  899. }
  900.  
  901. /*
  902.  * Interpolate given point list (defined via p_cntr) using Complete
  903.  * Spline interpolation.
  904.  */
  905. static complete_spline_interp(p_cntr, n, t_min, t_max, tx1, ty1, tx2, ty2)
  906. struct cntr_struct *p_cntr;
  907. int n;
  908. double t_min, t_max, tx1, ty1, tx2, ty2;
  909. {
  910.     double dt, *tangents_x, *tangents_y;
  911.     int i;
  912.  
  913.     tangents_x = (double *) alloc((unsigned) (sizeof(double) * n),
  914.                         "contour c_s_intr");
  915.     tangents_y = (double *) alloc((unsigned) (sizeof(double) * n),
  916.                         "contour c_s_intr");
  917.  
  918.     if (n > 1) prepare_spline_interp(tangents_x, tangents_y, p_cntr, n,
  919.                     t_min, t_max, tx1, ty1, tx2, ty2);
  920.     else {
  921.     free((char *) tangents_x);
  922.     free((char *) tangents_y);
  923.     return;
  924.     }
  925.  
  926.     dt = (t_max-t_min)/(n-1);
  927.  
  928.     add_cntr_point(p_cntr -> X, p_cntr -> Y);           /* First point. */
  929.  
  930.     for (i=0; i<n-1; i++) {
  931.         hermit_interp(p_cntr -> X, p_cntr -> Y,
  932.                      tangents_x[i], tangents_y[i],
  933.              p_cntr -> next -> X, p_cntr -> next -> Y,
  934.                      tangents_x[i+1], tangents_y[i+1], dt);
  935.  
  936.         p_cntr = p_cntr -> next;
  937.     }
  938.  
  939.     free((char *) tangents_x);
  940.     free((char *) tangents_y);
  941. }
  942.  
  943. /*
  944.  * Routine to calculate intermidiate value of the Hermit Blending function:
  945.  * This routine should be called only ONCE at the beginning of the program.
  946.  */
  947. static calc_hermit_table()
  948. {
  949.     int i;
  950.     double t, dt;
  951.  
  952.     hermit_table = (table_entry *) alloc ((unsigned) (sizeof(table_entry) *
  953.                         (num_approx_pts + 1)),
  954.                         "contour hermit table");
  955.     t = 0;
  956.     dt = 1.0/num_approx_pts;
  957.     for (i=0; i<=num_approx_pts; i++) {
  958.         hermit_table[i][0] = (t-1)*(t-1)*(2*t+1);             /* h00. */
  959.         hermit_table[i][1] = t*t*(-2*t+3);                 /* h10. */
  960.         hermit_table[i][2] = t*(t-1)*(t-1);                 /* h01. */
  961.         hermit_table[i][3] = t*t*(t-1);                     /* h11. */
  962.         t = t + dt;
  963.     }
  964. }
  965.  
  966. /*
  967.  * Routine to generate an hermit interpolation between two points given as
  968.  * two InterpStruct structures. Assume hermit_table is already calculated.
  969.  * Currently the points generated are printed to stdout as two reals (X, Y).
  970.  */
  971. static hermit_interp(x1, y1, tx1, ty1, x2, y2, tx2, ty2, dt)
  972. double x1, y1, tx1, ty1, x2, y2, tx2, ty2, dt;
  973. {
  974.     int i;
  975.     double x, y, vec_size, tang_size;
  976.  
  977.     tx1 *= dt;  ty1 *= dt; /* Normalize the tangents according to param. t.  */
  978.     tx2 *= dt;  ty2 *= dt;
  979.  
  980.     /* Normalize the tangents so that their magnitude will be 1/3 of the     */
  981.     /* segment length. This tumb rule guaranteed no cusps or loops!          */
  982.     /* Note that this normalization keeps continuity to be G1 (but not C1).  */
  983.     vec_size = sqrt(sqr(x1 - x2) + sqr(y2 - y1));
  984.     tang_size = sqrt(sqr(tx1) + sqr(ty1));                  /* Normalize T1. */
  985.     if (tang_size * 3 > vec_size) {
  986.     tx1 *= vec_size / (tang_size * 3);
  987.     ty1 *= vec_size / (tang_size * 3);
  988.     }
  989.     tang_size = sqrt(sqr(tx2) + sqr(ty2));                  /* Normalize T2. */
  990.     if (tang_size * 3 > vec_size) {
  991.     tx2 *= vec_size / (tang_size * 3);
  992.     ty2 *= vec_size / (tang_size * 3);
  993.     }
  994.  
  995.     for (i=1; i<=num_approx_pts; i++) {      /* Note we start from 1 - first */
  996.         x = hermit_table[i][0] * x1 +       /* point is not printed as it is */
  997.             hermit_table[i][1] * x2 +   /* redundent (last on last section). */
  998.             hermit_table[i][2] * tx1 +
  999.             hermit_table[i][3] * tx2;
  1000.         y = hermit_table[i][0] * y1 +
  1001.             hermit_table[i][1] * y2 +
  1002.             hermit_table[i][2] * ty1 +
  1003.             hermit_table[i][3] * ty2;
  1004.     add_cntr_point(x, y);
  1005.     }
  1006. }
  1007.  
  1008. /*
  1009.  * Routine to Set up the 3*N mat for solve_tri_diag routine used in the
  1010.  * Complete Spline Interpolation. Returns TRUE of calc O.K.
  1011.  * Gets the points list in p_cntr (Of length n) and with tangent vectors tx1,
  1012.  * ty1 at starting point and tx2, ty2 and end point.
  1013.  */
  1014. static prepare_spline_interp(tangents_x, tangents_y, p_cntr, n, t_min, t_max,
  1015.                tx1, ty1, tx2, ty2)
  1016. double tangents_x[], tangents_y[];
  1017. struct cntr_struct *p_cntr;
  1018. int n;
  1019. double t_min, t_max, tx1, ty1, tx2, ty2;
  1020. {
  1021.     int i;
  1022.     double *r, t, dt;
  1023.     tri_diag *m;           /* The tri-diagonal matrix is saved here. */
  1024.     struct cntr_struct *p;
  1025.  
  1026.     m = (tri_diag *) alloc((unsigned) (sizeof(tri_diag) * n),
  1027.                         "contour tri_diag");
  1028.     r = (double *) alloc((unsigned) (sizeof(double) * n),
  1029.                         "contour tri_diag2");
  1030.     n--;
  1031.  
  1032.     p = p_cntr;
  1033.     m[0][0] = 0.0;    m[0][1] = 1.0;    m[0][2] = 0.0;
  1034.     m[n][0] = 0.0;    m[n][1] = 1.0;    m[n][2] = 0.0;
  1035.     r[0] = tx1;                           /* Set start tangent. */
  1036.     r[n] = tx2;                         /* Set end tangent. */
  1037.     t = t_min;
  1038.     dt = (t_max-t_min)/n;
  1039.     for (i=1; i<n; i++) {
  1040.        t = t + dt;
  1041.        m[i][0] = dt;
  1042.        m[i][2] = dt;
  1043.        m[i][1] = 2 * (m[i][0] + m[i][2]);
  1044.        r[i] = m[i][0] * ((p -> next -> X) - (p -> X)) / m[i][2]
  1045.             + m[i][2] * ((p -> next -> next -> X) - 
  1046.                          (p -> next -> X)) / m[i][0];
  1047.        r[i] *= 3.0;
  1048.        p = p -> next;
  1049.     }
  1050.  
  1051.     if (!solve_tri_diag(m, r, tangents_x, n+1)) { /* Find the X(t) tangents. */
  1052.         free((char *) m);
  1053.         free((char *) r);
  1054.     int_error("Cannt interpolate X using complete splines", NO_CARET);
  1055.     }
  1056.  
  1057.     p = p_cntr;
  1058.     m[0][0] = 0.0;    m[0][1] = 1.0;    m[0][2] = 0.0;
  1059.     m[n][0] = 0.0;    m[n][1] = 1.0;    m[n][2] = 0.0;
  1060.     r[0] = ty1;                           /* Set start tangent. */
  1061.     r[n] = ty2;                         /* Set end tangent. */
  1062.     t = t_min;
  1063.     dt = (t_max-t_min)/n;
  1064.     for (i=1; i<n; i++) {
  1065.        t = t + dt;
  1066.        m[i][0] = dt;
  1067.        m[i][2] = dt;
  1068.        m[i][1] = 2 * (m[i][0] + m[i][2]);
  1069.        r[i] = m[i][0] * ((p -> next -> Y) - (p -> Y)) / m[i][2]
  1070.             + m[i][2] * ((p -> next -> next -> Y) -
  1071.                          (p -> next -> Y)) / m[i][0];
  1072.        r[i] *= 3.0;
  1073.        p = p -> next;
  1074.     }
  1075.  
  1076.     if (!solve_tri_diag(m, r, tangents_y, n+1)) { /* Find the Y(t) tangents. */
  1077.         free((char *) m);
  1078.         free((char *) r);
  1079.     int_error("Cannt interpolate Y using complete splines", NO_CARET);
  1080.     }
  1081.     free((char *) m);
  1082.     free((char *) r);
  1083. }
  1084.  
  1085. /*
  1086.  * Solve tri diagonal linear system equation. The tri diagonal matrix is
  1087.  * defined via matrix M, right side is r, and solution X i.e. M * X = R.
  1088.  * Size of system given in n. Return TRUE if solution exist.
  1089.  */
  1090. static int solve_tri_diag(m, r, x, n)
  1091. tri_diag m[];
  1092. double r[], x[];
  1093. int n;
  1094. {
  1095.     int i;
  1096.     double t;
  1097.  
  1098.     for (i=1; i<n; i++) {   /* Eliminate element m[i][i-1] (lower diagonal). */
  1099.     if (m[i-1][1] == 0) return FALSE;
  1100.     t = m[i][0] / m[i-1][1];        /* Find ratio between the two lines. */
  1101.     m[i][0] = m[i][0] - m[i-1][1] * t;
  1102.     m[i][1] = m[i][1] - m[i-1][2] * t;
  1103.     r[i] = r[i] - r[i-1] * t;
  1104.     }
  1105.     /* Now do back subtitution - update the solution vector X: */
  1106.     if (m[n-1][1] == 0) return FALSE;
  1107.     x[n-1] = r[n-1] / m[n-1][1];               /* Find last element. */
  1108.     for (i=n-2; i>=0; i--) {
  1109.     if (m[i][1] == 0) return FALSE;
  1110.     x[i] = (r[i] - x[i+1] * m[i][2]) / m[i][1];
  1111.     }
  1112.     return TRUE;
  1113. }
  1114.  
  1115. /*
  1116.  * Generate a Bspline curve defined by all the points given in linked list p:
  1117.  * Algorithm: using deBoor algorithm
  1118.  * Note: if Curvekind is OPEN_CONTOUR than Open end knot vector is assumed,
  1119.  *       else (CLOSED_CONTOUR) Float end knot vector is assumed.
  1120.  * It is assumed that num_of_points is at list 2, and order of Bspline is less
  1121.  * than num_of_points!
  1122.  */
  1123. static gen_bspline_approx(p_cntr, num_of_points, order, contour_kind)
  1124. struct cntr_struct *p_cntr;
  1125. int num_of_points, order, contour_kind;
  1126. {
  1127.     int i, knot_index = 0, pts_count = 1;
  1128.     double dt, t, next_t, t_min, t_max, x, y;
  1129.     struct cntr_struct *pc_temp = p_cntr, *pc_tail;
  1130.  
  1131.     /* If the contour is Closed one we must update few things:               */
  1132.     /* 1. Make the list temporary circular, so we can close the contour.     */
  1133.     /* 2. Update num_of_points - increase it by "order-1" so contour will be */
  1134.     /*    closed. This will evaluate order more sections to close it!        */
  1135.     if (contour_kind == CLOSED_CONTOUR) {
  1136.     pc_tail = p_cntr;
  1137.     while (pc_tail -> next) pc_tail = pc_tail -> next;/* Find last point.*/
  1138.     pc_tail -> next = p_cntr;   /* Close contour list - make it circular.*/
  1139.     num_of_points += order;
  1140.     }
  1141.  
  1142.     /* Find first (t_min) and last (t_max) t value to eval: */
  1143.     t = t_min = fetch_knot(contour_kind, num_of_points, order, order);
  1144.     t_max = fetch_knot(contour_kind, num_of_points, order, num_of_points);
  1145.     next_t = t_min + 1.0;
  1146.     knot_index = order;
  1147.     dt = 1.0/num_approx_pts;            /* Number of points per one section. */
  1148.  
  1149.  
  1150.     while (t<t_max) {
  1151.     if (t > next_t) {
  1152.         pc_temp = pc_temp -> next;     /* Next order ctrl. pt. to blend. */
  1153.             knot_index++;
  1154.         next_t += 1.0;
  1155.     }
  1156.         eval_bspline(t, pc_temp, num_of_points, order, knot_index,
  1157.                         contour_kind, &x, &y);   /* Next pt. */
  1158.     add_cntr_point(x, y);
  1159.     pts_count++;
  1160.     /* As we might have some real number round off problems we must      */
  1161.     /* test if we dont produce too many points here...                   */
  1162.     if (pts_count + 1 == num_approx_pts * (num_of_points - order) + 1)
  1163.             break;
  1164.         t += dt;
  1165.     }
  1166.  
  1167.     eval_bspline(t_max - EPSILON, pc_temp, num_of_points, order, knot_index,
  1168.         contour_kind, &x, &y);
  1169.     /* If from round off errors we need more than one last point: */
  1170.     for (i=pts_count; i<num_approx_pts * (num_of_points - order) + 1; i++)
  1171.     add_cntr_point(x, y);                /* Complete the contour. */
  1172.  
  1173.     if (contour_kind == CLOSED_CONTOUR)     /* Update list - un-circular it. */
  1174.     pc_tail -> next = NULL;
  1175. }
  1176.  
  1177. /*
  1178.  * The recursive routine to evaluate the B-spline value at point t using
  1179.  * knot vector PKList, and the control points Pdtemp. Returns x, y after the
  1180.  * division by the weight w. Note that Pdtemp points on the first control
  1181.  * point to blend with. The B-spline is of order order.
  1182.  */
  1183. static eval_bspline(t, p_cntr, num_of_points, order, j, contour_kind, x, y)
  1184. double t;
  1185. struct cntr_struct *p_cntr;
  1186. int num_of_points, order, j, contour_kind;
  1187. double *x, *y;
  1188. {
  1189.     int i, p;
  1190.     double ti, tikp, *dx, *dy;      /* Copy p_cntr into it to make it faster. */
  1191.  
  1192.     dx = (double *) alloc((unsigned) (sizeof(double) * (order + j)),
  1193.                         "contour b_spline");
  1194.     dy = (double *) alloc((unsigned) (sizeof(double) * (order + j)),
  1195.                         "contour b_spline");
  1196.     /* Set the dx/dy - [0] iteration step, control points (p==0 iterat.): */
  1197.     for (i=j-order; i<=j; i++) {
  1198.         dx[i] = p_cntr -> X;
  1199.         dy[i] = p_cntr -> Y;
  1200.         p_cntr = p_cntr -> next;
  1201.     }
  1202.  
  1203.     for (p=1; p<=order; p++) {        /* Iteration (b-spline level) counter. */
  1204.     for (i=j; i>=j-order+p; i--) {           /* Control points indexing. */
  1205.             ti = fetch_knot(contour_kind, num_of_points, order, i);
  1206.             tikp = fetch_knot(contour_kind, num_of_points, order, i+order+1-p);
  1207.         if (ti == tikp) {   /* Should not be a problems but how knows... */
  1208.         }
  1209.         else {
  1210.         dx[i] = dx[i] * (t - ti)/(tikp-ti) +         /* Calculate x. */
  1211.             dx[i-1] * (tikp-t)/(tikp-ti);
  1212.         dy[i] = dy[i] * (t - ti)/(tikp-ti) +         /* Calculate y. */
  1213.             dy[i-1] * (tikp-t)/(tikp-ti);
  1214.         }
  1215.     }
  1216.     }
  1217.     *x = dx[j]; *y = dy[j];
  1218.     free((char *) dx);
  1219.     free((char *) dy);
  1220. }
  1221.  
  1222. /*
  1223.  * Routine to get the i knot from uniform knot vector. The knot vector
  1224.  * might be float (Knot(i) = i) or open (where the first and last "order"
  1225.  * knots are equal). contour_kind determines knot kind - OPEN_CONTOUR means
  1226.  * open knot vector, and CLOSED_CONTOUR selects float knot vector.
  1227.  * Note the knot vector is not exist and this routine simulates it existance
  1228.  * Also note the indexes for the knot vector starts from 0.
  1229.  */
  1230. static double fetch_knot(contour_kind, num_of_points, order, i)
  1231. int contour_kind, num_of_points, order, i;
  1232. {
  1233.     switch (contour_kind) {
  1234.     case OPEN_CONTOUR:
  1235.         if (i <= order) return 0.0;
  1236.         else if (i <= num_of_points) return (double) (i - order);
  1237.          else return (double) (num_of_points - order);
  1238.     case CLOSED_CONTOUR:
  1239.         return (double) i;
  1240.     default: /* Should never happen */
  1241.         return 1.0;
  1242.     }
  1243. }
  1244.